Sesión 4

Curso: R para usuarios de Stata

Percy Soto-Becerra, M.D., M.Sc(c)

InkaStats Academy | Medical Branch

2023

https://github.com/psotob91

Inferencia estadística

Inferencia estadística

Inferencia en Stata vs. R

  • Stata tiene varios comandos oficiales para realizar pruebas de hipótesis o estimar intervalos de confianza de manera clásica.

  • Tiene comandos para pruebas:

    • Paramétricas
    • No paramétricas basadas en rangos
    • No paramétricas basadas en remuestreo

  • R base también tiene funciones para realizar inferencia estadística clásica.

  • Problema: No consistencia y dificulta de reporte reproducible.

  • Otros paquetes suplen este problema:

    • {rstatix}
    • {infer}

Paquete rstatix

  • Provee marco simple e intuitivo compatible con tidyverse y pipelines.

  • Permite realizar varios tipos de pruebas:

    • Pruebas t, prueba de Wilcoxon, ANOVA, Kruskal-Wallis y análisis de correlación.

  • Ingresa data.frame y retorna un data.frame

  • Múltiples funciones adicionales para remodelar, reordenar, manipular y visualizar los resultados de las pruebas.

  • Descargar e instalar paquete:

#install.packages(rstatix)
library(rstatix)

Prueba t de Student para un grupo

ttest y == numero

t_test(y ~ 1, mu = numero, detailed = TRUE)

Ejemplo

  • Se usa función t_test() de paquete {rstatix}.

  • Permite obtener el valor p para una hipótesis dada y el intervalo de confianza para la media.

datos %>% 
  t_test(weight ~ 1, mu = 100, detailed = TRUE) 
# A tibble: 1 × 12
  estimate .y.    group1 group2         n stati…¹        p    df conf.…² conf.…³
*    <dbl> <chr>  <chr>  <chr>      <int>   <dbl>    <dbl> <dbl>   <dbl>   <dbl>
1     60.5 weight 1      null model    50   -35.6 1.15e-36    49    58.3    62.8
# … with 2 more variables: method <chr>, alternative <chr>, and abbreviated
#   variable names ¹​statistic, ²​conf.low, ³​conf.high

  • El argumento mu = número es un valor de referencia contra el que deseamos compararnos.

  • El argumento detailed = TRUE permite obtener también los intervalos de confianza.

  • La función gt() es para mejorar visualización del resultado.

library(gt)
datos %>% 
  t_test(weight ~ 1, mu = 100, detailed = TRUE) %>% 
  gt()
estimate .y. group1 group2 n statistic p df conf.low conf.high method alternative
60.524 weight 1 null model 50 -35.57465 1.15e-36 49 58.29404 62.75396 T-test two.sided

  • Agregaremos gt() a partir de ahora.

  • Interpretación:

La media estimada de glucosa fue de 95.48 mg/dL (IC95% 71.6 mg/dL a 119.4 mg/dL). Aunque este valor fue menor a 100 mg/dL (valor de referencia), no es posible concluir que la media poblacional es diferente del valor de refrencia 100 mg/dL (p = 0.627).

¿El nivel medio de estradiol es diferente de 50 UI en la población?

Para variable estradiol:

datos %>% 
  t_test(e2 ~ 1, mu = 50, detailed = TRUE) %>% 
  gt()
estimate .y. group1 group2 n statistic p df conf.low conf.high method alternative
112.6717 e2 1 null model 53 6.459575 3.6e-08 52 93.20293 132.1405 T-test two.sided

Prueba del signo para un grupo

*No equivalente oficial aún.

sign_test(y ~ 1, mu = numero, detailed = TRUE)

Pruebas de hipotesis en Stata y rstatix de R

La estructura es nombre_de_prueba_test

  • ttest …, paired
  • signtest …
  • signrank …
  • kwallis …
  • oneway …
  • anova …

  • t_test(…)
  • sign_test(…)
  • wilcox_test(…)
  • kruskal_test(…)
  • anova_test(…)

Ejemplo de prueba del signo para un grupo

Prueba del Signo:

  • Permite comparar la mediana (ya no la media) contra un valor de referencia.

  • No permite estimar intervalos de confianza.

datos %>% 
  sign_test(weight ~ 1, mu = 0, detailed = TRUE) %>% 
  gt()
estimate .y. group1 group2 n statistic p df conf.low conf.high method alternative
58.8 weight 1 null model 50 50 1.78e-15 50 56 62.5 Sign-Test two.sided
  • Interpretación:

La mediana estimada de glucosa fue de 95.53 mg/dL. Aunque este valor fue menor a 100 mg/dL (valor de referencia), no es posible concluir que la media poblacional es diferente del valor de refrencia (p = 0.627).

¿El nivel mediano de estradiol es diferente de 50 UI en la población?

datos %>% 
  sign_test(e2 ~ 1, mu = 50, detailed = TRUE) %>% 
  gt()
estimate .y. group1 group2 n statistic p df conf.low conf.high method alternative
103.95 e2 1 null model 53 42 2.25e-05 53 75.78 122.97 Sign-Test two.sided
  • Interpretación:

La mediana estimada de estradiol fue de 103.95 UI. Este valor fue mayor a 50 UI (valor de referencia) (p < 0.001).

Estimación y PH para dos grupos dependientes

  • Se usa la prueba t de Student para datos pareados.

  • Como los datos están pareados 1:1, la prueba trabaja sobre la diferencia de los valores individuales convirtiendo el problema en una sola muestra.

  • Por lo tanto, los supuestos son los mismos que para la `prueba t de Student de un solo grupo.

Ejemplo: ¿Cuál es el efecto de la Vitamina C en el crecimiento de los dientes en cuyes experimentales?

len: Longitud de dientes al final del experimento.

supp: Suplemento de vitamina C (VC vs. OJ)

Error in file(file, "rt"): no se puede abrir la conexión
Error in mutate(., dif = VC - OJ): objeto 'df2' no encontrado
df %>% t_test(len ~ supp, 
              paired = TRUE, 
              detailed = TRUE) %>% gt()
estimate .y. group1 group2 n1 n2 statistic p df conf.low conf.high method alternative
3.7 len OJ VC 30 30 3.302585 0.00255 29 1.408659 5.991341 T-test two.sided

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Normalidad o cumplimiento del TLC.

library(ggpubr)
df3 %>% 
  ggqqplot(x = "dif")
Error in ggqqplot(., x = "dif"): objeto 'df3' no encontrado

  • Variable al menos en escala ordinal

  • Supuestos implícitos (no error de medición ni sesgos)

  • Se tienen dos opciones: Prueba del signo vs. Prueba de rangos signados de Wilcoxon.

  • Se trabaja igual que para t-test pareada, con al diferencia de los valores de tal forma que se tiene solo un grupo.

  • Prueba del signo no requiere supuesto de simetría de distribución para comparar medianas.

  • Prueba de Wilcoxon requiere supuesto de simetría de distribución para comparar medianas.

Prueba del signo:

df %>% sign_test(len ~ supp, 
              detailed = TRUE) %>% gt()
estimate .y. group1 group2 n1 n2 statistic p df conf.low conf.high method alternative
4.85 len OJ VC 30 30 19 0.136 29 -0.3 7.9 Sign-Test two.sided

Prueba de Wilcoxon:

df %>% wilcox_test(len ~ supp, 
              detailed = TRUE) %>% gt()
estimate .y. group1 group2 n1 n2 statistic p conf.low conf.high method alternative
4.000014 len OJ VC 30 30 575.5 0.0645 -0.1000097 8.500017 Wilcoxon two.sided

Supuestos de la Prueba del Signo

  • Aleatorización

  • Indenpendencia de observaciones

  • Escala de medición al menos ordinal

Supuestos de la Prueba de Wilcoxon

  • Aleatorización

  • Indenpendencia de observaciones

  • Escala de medición al menos ordinal

  • Distribución simétrica

Estimación y PH para dos grupos independientes

  • Se usa prueba t de Student para grupos independientes.
datos2 %>% 
  t_test(weight ~ treat, detailed = TRUE, var.equal = TRUE) %>% 
  gt()
estimate estimate1 estimate2 .y. group1 group2 n1 n2 statistic p df conf.low conf.high method alternative p.adj p.adj.signif
-2.2812500 59.21875 61.50000 weight Placebo Dosis 1 16 16 -0.7745439 0.445 30 -8.296318 3.733818 T-test two.sided 1 ns
-1.5979167 59.21875 60.81667 weight Placebo Dosis 2 16 18 -0.7126840 0.481 32 -6.164947 2.969114 T-test two.sided 1 ns
0.6833333 61.50000 60.81667 weight Dosis 1 Dosis 2 16 18 0.2249600 0.823 32 -5.504008 6.870674 T-test two.sided 1 ns

Supuestos de la Prueba t de Student para grupos independientes

  • Aleatorización

  • Indenpendencia de observaciones

  • Normalidad

  • Homogeneidad de varianzas (*)

datos2 %>% 
  ggqqplot("weight", color = "treat")

datos2 %>% 
  group_by(treat) %>% 
  skim(weight)
Data summary
Name Piped data
Number of rows 53
Number of columns 14
_______________________
Column type frequency:
numeric 1
________________________
Group variables treat

Variable type: numeric

skim_variable treat n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
weight Placebo 1 0.94 59.22 5.61 50.1 54.90 58.55 61.12 71.0 ▅▇▆▃▃
weight Dosis 1 1 0.94 61.50 10.36 50.1 55.68 59.30 64.22 92.1 ▇▇▁▁▁
weight Dosis 2 1 0.95 60.82 7.24 51.4 54.50 60.05 65.88 74.5 ▇▂▃▃▃
  • Si no confiamos en la homogeneidad de varianzas, podemos hacer la prueba robusta a la heterogeneidad de varianzas usando la corrección de Satterwhite:
datos2 %>% 
  t_test(weight ~ treat, detailed = TRUE, var.equal = FALSE) %>% 
  gt()
estimate estimate1 estimate2 .y. group1 group2 n1 n2 statistic p df conf.low conf.high method alternative p.adj p.adj.signif
-2.2812500 59.21875 61.50000 weight Placebo Dosis 1 16 16 -0.7745439 0.446 23.09487 -8.372645 3.810145 T-test two.sided 1 ns
-1.5979167 59.21875 60.81667 weight Placebo Dosis 2 16 18 -0.7235647 0.475 31.45480 -6.099330 2.903497 T-test two.sided 1 ns
0.6833333 61.50000 60.81667 weight Dosis 1 Dosis 2 16 18 0.2203118 0.827 26.44841 -5.686973 7.053640 T-test two.sided 1 ns
  • Se usa prueba de suma de rangos de Wilcoxon, también denominada U de Mann Withney.
datos2 %>% 
  wilcox_test(e2 ~ treat, detailed = TRUE) %>% 
  gt()
estimate .y. group1 group2 n1 n2 statistic p conf.low conf.high method alternative p.adj p.adj.signif
-32.14 e2 Placebo Dosis 1 17 17 106 0.193 -64.80 19.98 Wilcoxon two.sided 0.579 ns
-20.27 e2 Placebo Dosis 2 17 19 125 0.257 -68.78 26.52 Wilcoxon two.sided 0.579 ns
-0.40 e2 Dosis 1 Dosis 2 17 19 160 0.975 -41.24 33.54 Wilcoxon two.sided 0.975 ns

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Escala de medición al menos ordinal

  • Igual forma de distribución entre los grupos (solo diferencia en la posición)

Estimación y Pruebas de Hipótesis para Variables Respuesta Numéricas: Tres o más grupos

ANOVA one-way en R

library(rstatix)
datos %>% 
  anova_test(weight ~ treat, detailed = TRUE)
ANOVA Table (type II tests)

  Effect    SSn      SSd DFn DFd     F     p p<.05   ges
1  treat 49.021 7171.016   2 100 0.342 0.711       0.007
  • La prueba ANOVA one-way no genera resultados puntuales, debemos generarlos con otra función:
datos %>% 
  group_by(treat) %>% 
  get_summary_stats(weight, type = "mean_sd")
# A tibble: 3 × 5
  treat   variable     n  mean    sd
  <chr>   <fct>    <dbl> <dbl> <dbl>
1 Dosis 1 weight      33  62.1 10.4 
2 Dosis 2 weight      37  61.2  8.13
3 Placebo weight      33  60.4  6.40
  • También es bueno visualizar los patrones mediante gráfico de cajas:
library(ggpubr)
datos %>% 
  ggboxplot(x = "treat", y = "weight")

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Normalidad

    • Supuesto de normalidad condicional: \(Y_i \sim Normal\), es decir \(Y|x \sim Normal\)
library(ggpubr)
datos %>% 
  ggqqplot("weight", facet.by = "treat")

+ Una mejor forma de evaluar la normalidad es usando los residuos.

  • Homogeneidad de varianzas (*)

    • Se podría evaluar si las rectas de los gráficos qqplot normal son paralelas o aproximadamente paralelas:
datos %>% 
  ggqqplot("weight", color = "treat")

  • Otra opción es evaluar los residuales. Esto lo veremos en regresión lineal.

  • Es preferible asumir que no hay homocedaticidad y realizar un ajuste robusto a heterogeneidad de varianzas.

datos %>% 
  anova_test(weight ~ treat, detailed = TRUE, white.adjust = TRUE)
ANOVA Table (type II tests)

  Effect DFn DFd     F    p p<.05
1  treat   2 100 0.329 0.72
  • No valores extremos.

(*) Este supuesto puede flexibilizarse:

  1. Usando algún método robusto de estimación de varianza.

  2. Usando remuestreo (Prueba de permutación o Bootstrapping)

datos %>% 
  anova_test(lh ~ treat, detailed = TRUE)
ANOVA Table (type II tests)

  Effect     SSn      SSd DFn DFd     F    p p<.05   ges
1  treat 470.911 21640.35   2 103 1.121 0.33       0.021
  • La prueba ANOVA one-way no genera resultados puntuales, debemos generarlos con otra función:
datos %>% 
  group_by(treat) %>% 
  get_summary_stats(lh, type = "mean_sd")
# A tibble: 3 × 5
  treat   variable     n  mean    sd
  <chr>   <fct>    <dbl> <dbl> <dbl>
1 Dosis 1 lh          34 13.7  16.6 
2 Dosis 2 lh          38  8.65  7.75
3 Placebo lh          34 11.6  17.7
  • También es bueno visualizar los patrones mediante gráfico de cajas:
datos %>% 
  ggboxplot(x = "treat", y = "lh")

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Normalidad

    • Supuesto de normalidad condicional: \(Y_i \sim Normal\), es decir \(Y|x \sim Normal\)
library(ggpubr)
datos %>% 
  ggqqplot("lh", facet.by = "treat")

+ Una mejor forma de evaluar la normalidad es usando los residuos.

  • Homogeneidad de varianzas (*)

    • Se podría evaluar si las rectas de los gráficos qqplot normal son paralelas o aproximadamente paralelas:
datos %>% 
  ggqqplot("lh", color = "treat")

  • Otra opción es evaluar los residuales. Esto lo veremos en regresión lineal.

  • A veces es preferible asumir que no hay homocedaticidad y realizar un ajuste robusto a heterogeneidad de varianzas.

datos %>% 
  anova_test(lh ~ treat, detailed = TRUE, white.adjust = TRUE)
ANOVA Table (type II tests)

  Effect DFn DFd     F     p p<.05
1  treat   2 103 1.502 0.227
  • No valores extremos.

(*) Este supuesto puede flexibilizarse:

  1. Usando algún método robusto de estimación de varianza.

  2. Usando remuestreo (Prueba de permutación o Bootstrapping)

Kruskal Wallis test

datos %>% 
  kruskal_test(weight ~ treat)
# A tibble: 1 × 6
  .y.        n statistic    df     p method        
* <chr>  <int>     <dbl> <int> <dbl> <chr>         
1 weight   106    0.0347     2 0.983 Kruskal-Wallis
  • La prueba Kruskal Wallis no genera resultados puntuales, debemos generarlos con otra función:
datos %>% 
  group_by(treat) %>% 
  get_summary_stats(weight, type = "median_iqr")
# A tibble: 3 × 5
  treat   variable     n median   iqr
  <chr>   <fct>    <dbl>  <dbl> <dbl>
1 Dosis 1 weight      33     60   9  
2 Dosis 2 weight      37     61  12  
3 Placebo weight      33     59   8.4
  • También es bueno visualizar las distribuciones
library(ggpubr)
datos %>% 
  ggboxplot(y = "weight", x = "treat")

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Misma distribución salvo por la mediana:

    • Esto significa que Kruskal Wallis también es perjudicado en cierto modo por la heterogeneidad de varianzas porque afecta el supuesto de igualdad de distribuciones

    • Podemos usar gráfico de densidades:

library(ggpubr)
datos %>% 
  ggdensity(x = "weight", 
            y = "..density..", 
            fill = "treat", 
            color = "treat", 
            add = "median")

+ Podemos usar gráfico de violin: Combina cajas y densidad
library(ggpubr)
datos %>% 
  ggviolin(x = "treat", 
            y = "weight", 
            color = "treat", 
            add = "boxplot")

Comparaciones múltiples

  • Hay varios enfoques de comparaciones múltiples.

  • Algunas son pre-hoc y otras post-hoc.

  • Veremos solo un par para ejemplificar:

  • Para ANOVA one-way clásico:
aov(weight ~ treat, data = datos) %>% 
  tukey_hsd()
# A tibble: 3 × 9
  term  group1  group2  null.value estimate conf.low conf.high p.adj p.adj.sig…¹
* <chr> <chr>   <chr>        <dbl>    <dbl>    <dbl>     <dbl> <dbl> <chr>      
1 treat Dosis 1 Dosis 2          0   -0.974    -5.80      3.85 0.881 ns         
2 treat Dosis 1 Placebo          0   -1.72     -6.68      3.24 0.689 ns         
3 treat Dosis 2 Placebo          0   -0.745    -5.57      4.08 0.928 ns         
# … with abbreviated variable name ¹​p.adj.signif
  • Si homogeneidad de varianzas no se cumple:
datos %>% 
  games_howell_test(weight ~ treat)
# A tibble: 3 × 8
  .y.    group1  group2  estimate conf.low conf.high p.adj p.adj.signif
* <chr>  <chr>   <chr>      <dbl>    <dbl>     <dbl> <dbl> <chr>       
1 weight Dosis 1 Dosis 2   -0.974    -6.39      4.44 0.903 ns          
2 weight Dosis 1 Placebo   -1.72     -6.85      3.42 0.701 ns          
3 weight Dosis 2 Placebo   -0.745    -4.91      3.43 0.904 ns
  • En epidemiología e investigación clínica se prefiere ajustes más conservadores. Uno de ellos es el ajuste de Bonferroni (uno de los más conservadores)
library(emmeans)
datos %>% 
  emmeans_test(weight ~ treat, p.adjust.method = "bonferroni")
# A tibble: 3 × 9
  term  .y.    group1  group2     df statistic     p p.adj p.adj.signif
* <chr> <chr>  <chr>   <chr>   <dbl>     <dbl> <dbl> <dbl> <chr>       
1 treat weight Dosis 1 Dosis 2   100     0.480 0.632     1 ns          
2 treat weight Dosis 1 Placebo   100     0.824 0.412     1 ns          
3 treat weight Dosis 2 Placebo   100     0.367 0.714     1 ns          
datos %>% 
  dunn_test(weight ~ treat)
# A tibble: 3 × 9
  .y.    group1  group2     n1    n2 statistic     p p.adj p.adj.signif
* <chr>  <chr>   <chr>   <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 weight Dosis 1 Dosis 2    33    37   -0.0300 0.976     1 ns          
2 weight Dosis 1 Placebo    33    33   -0.173  0.863     1 ns          
3 weight Dosis 2 Placebo    37    33   -0.148  0.882     1 ns

Pruebas de hipótesis para variables respuesta categóricas

Prueba Chi cuadrado

table(datos2$married2, datos2$treat) -> xtab

xtab
                
                 Placebo Dosis 1 Dosis 2
  Without couple       9      10       8
  With couple          8       7      11
xtab  %>% 
  chisq_test()
# A tibble: 1 × 6
      n statistic     p    df method          p.signif
* <int>     <dbl> <dbl> <int> <chr>           <chr>   
1    53      1.04 0.594     2 Chi-square test ns      
xtab  %>% 
  chisq_test() %>%  
  chisq_descriptives()
# A tibble: 6 × 9
  Var1           Var2    observed  prop row.prop col.prop expec…¹  resid std.r…²
  <fct>          <fct>      <int> <dbl>    <dbl>    <dbl>   <dbl>  <dbl>   <dbl>
1 Without couple Placebo        9 0.170    0.333    0.529    8.66  0.115   0.200
2 With couple    Placebo        8 0.151    0.308    0.471    8.34 -0.118  -0.200
3 Without couple Dosis 1       10 0.189    0.370    0.588    8.66  0.455   0.789
4 With couple    Dosis 1        7 0.132    0.269    0.412    8.34 -0.464  -0.789
5 Without couple Dosis 2        8 0.151    0.296    0.421    9.68 -0.540  -0.962
6 With couple    Dosis 2       11 0.208    0.423    0.579    9.32  0.550   0.962
# … with abbreviated variable names ¹​expected, ²​std.resid
datos2 %>% 
  tbl_cross(
    row = married2,
    col = treat, 
    percent = "row"
  ) %>% 
  add_p()
treat Total p-value1
Placebo Dosis 1 Dosis 2
Marital status, recat 0.6
    Without couple 9 (33%) 10 (37%) 8 (30%) 27 (100%)
    With couple 8 (31%) 7 (27%) 11 (42%) 26 (100%)
Total 17 (32%) 17 (32%) 19 (36%) 53 (100%)
1 Pearson's Chi-squared test

Prueba exacta de Fisher

table(datos2$married2, datos2$treat) -> xtab

xtab
                
                 Placebo Dosis 1 Dosis 2
  Without couple       9      10       8
  With couple          8       7      11
xtab  %>% 
  fisher_test()
# A tibble: 1 × 3
      n     p p.signif
* <int> <dbl> <chr>   
1    53 0.696 ns      
datos2 %>% 
  tbl_cross(
    row = married2,
    col = treat, 
    percent = "row"
  ) %>% 
  add_p(
     test = "fisher.test"
  )
treat Total p-value1
Placebo Dosis 1 Dosis 2
Marital status, recat 0.7
    Without couple 9 (33%) 10 (37%) 8 (30%) 27 (100%)
    With couple 8 (31%) 7 (27%) 11 (42%) 26 (100%)
Total 17 (32%) 17 (32%) 19 (36%) 53 (100%)
1 Fisher's exact test

Creación de tablas reproducibles con gtsummary

  • Se puede agregar el valor p en las tablas de gtsummary:
library(gtsummary)
datos2 %>%
  dplyr::select(treat, age, weight, e2, fsh, lh, race, married) %>% 
  tbl_summary(
    by = treat
  ) %>%  
  add_p()
Characteristic Placebo, N = 171 Dosis 1, N = 171 Dosis 2, N = 191 p-value2
Age, years 34.0 (32.0, 37.0) 31.0 (27.0, 38.0) 35.0 (30.0, 36.5) >0.9
Weight, kg 59 (55, 61) 59 (56, 64) 60 (54, 66) >0.9
    Unknown 1 1 1
Estradiol 71 (40, 134) 112 (91, 130) 110 (76, 158) 0.4
Folicullo stimulant hormon 3.66 (2.17, 5.49) 4.00 (2.61, 4.83) 4.37 (2.54, 5.30) >0.9
Luteinizant Hormon 6 (3, 15) 7 (4, 15) 6 (4, 14) >0.9
Race
    Mestiza 17 (100%) 17 (100%) 19 (100%)
Marital status 0.6
    Single 8 (47%) 9 (53%) 8 (42%)
    Divorced 0 (0%) 1 (5.9%) 0 (0%)
    Widowed 1 (5.9%) 0 (0%) 0 (0%)
    Cohabiting 0 (0%) 1 (5.9%) 0 (0%)
    Married 8 (47%) 6 (35%) 11 (58%)
1 Median (IQR); n (%)
2 Kruskal-Wallis rank sum test; Fisher's exact test

Recomendaciones de reporte

  • No reporte valores p innecesarios, que no respondan sus preguntas de investigación pre-definidas.

  • Es una (muy difundida) mala práctica reportar tablas comparativas descriptivas con valores p.

  • Varias guías, comenzando por STROBE (para observacionales) y CONSORT (para ensayos clínicos), explícitamente recomiendan en contra de reportar estas tablas.

    • Solo haga una comparación descriptiva de los resultados de la tabla.

      • Por ejemplo, 14% vs. 20%, diferencia de 6% es mejor que “diferencia estadísticamente significativa”.

    • Reservese la inferencia para responder la pregunta principal:

      • Ahorra grados de libertad y previene grados de libertad "fantasma" (‘ghost degree of freedom’ según Frank Harrell)

¿Qué dice STROBE?

  • La guía “Strengthening the reporting of observational studies in epidemiology (STROBE) statement” da algunas recomendaciones para tablas tipo 1:

¿Qué dice STROBE?

  • STROBE recomienda en contra de reportar valores p en las tablas descriptivas!!

Fuente: Enlace aquí

  • Lamentablemente sigue siendo una mala práctica muy difundida: ¡Evitemosla en estudios observacionales!

Modelos de Regresión

Modelos de Regresión

Modelos de regresión

Simple

comando_modelo y x

funcion_modelo(y ~ x)

Múltiple

comando_modelo y x1 x2

funcion_modelo(y ~ x1 + x2, data = datos)

Notación fórmulas en regresión

Stata R
y varnum varcat y ~ varnum + varcat
y varnum, nocons y ~ 0 + varnum
y i.varcat y ~ varcat
y c.varnum1#c.varnum2 y ~ varnum1:varnum2
y c.varnum1##c.varnum2 y ~ varnum1*varnum2
y c.varnum1##i.varcat y ~ varnum1*varcat

Regresión lineal de Mínimos Cuadrados Ordinarios


regress y x1 x2 c.x3##i.x4 [weight = pesos]

modelo <- lm(y ~ x1 + x2 + x3*x4, data = datos, weights = pesos)

summary(modelo) # Opción de R base

modelo <- lm(y ~ x1 + x2 + x3*x4, data = datos, weights = pesos)

tidy(modelo) # Opción de R tidyverse con paquete broom

GLM

 


glm y x1 x2 [weight = pond], family(familia) link(enlace)

mod <- glm(y ~ x1 + x2,
            family = familia(link = “enlace”),
            data = datos, weights = pond)

summary(mod) # Opción de R base

tidy(mod) # Opción de R tidyverse con paquete broom

GLM family - link function

Stata Stata glm equivalente R
regress family(gaussian) link(identity) family = gaussian(link = “identity”)
logit family(binomial) link(logit) family = binomial(link = “logit”)
poisson family(poisson) link(log) family = poisson(link = “log”)
nbreg family(nbinomial ml) link(log) Otra función: glm.nb()
NA family(gamma) link(power -1) family = Gamma(link = “inverse”)
NA family(igaussian) link(power -2) family = inverse.gaussian(link = “1/mu^2”)
NA NA family = quasi(link = “identity”, variance = “constant”)
NA NA family = quasibinomial(link = “logit”)
NA NA family = quasipoisson(link = “log”)

Post-estimacion

Explorar resultados post-estimación

ereturn list

names(mod) # Aunque lo mejor es explorar la lista guardada

Generar residuos

predict cualq_nombre, tipo_residuo

cualq_nombre <- residuals(mod, type = “tipo_residuo”)

cualq_nombre <- augment(mod) # Genera lista de residuos y otros valores de interés

Supuestos

  • Paso 1: Crear los residuos y otros indicadores de interés.
  • Paso 2: Generar las gráficas con estos insumos.

  • Usar un wrapper.

    • Stata tiene solo un empaquetado para regresión lineal con regress, pero aún no tiene empaquetados oficiales para glm u otras regresiones:

      • Desarrollé un paquete inspirado en R, glm_diagnostic (Enlace aquí)
      • Aún en versión alfa con planes para extenderse a modelos no glm (cero inflados, beta y fraccionales, Cox y sobrevida paramétrica) y subirlo al repo oficial de Stata.

    • R tiene infinidad de empaquetados, los más populares:

      • plot() de R base
      • autoplot() de {ggfortify}
      • check_model() de {performance}
      • Otro basados en paquete {car}

Reporte reproducible

  • Se puede usar collect en compañía de putdocx, putpdf o putexcel.

  • También hay algunos empaquetados populares de la era previa a collect:

  • Se puede usar broom para obtener resultados en data.frames y pipearlos dentro de paquetes para crear tablas personalizadas: flextable, gt, huxtable, kable.

    • Pro: Tablas super-personalizadas / Contra: Mucho más código, similar a collect.

  • Se pueden usar empaquetados!

Empaquetados de R

  • Realmente hay infinidad de empaquetados (wrappers) en R para tablas reproducibles.

  • Los más populares, estables y confiables son:

    • finalfit
    • stargazer
    • tbl_uvregression y tbl_regression

¡Gracias!
¿Preguntas?




https://github.com/psotob91

percys1991@gmail.com

Curso: R para usuarios de Stata - Ed. 1 (2023)